quietly {
	
**********
*rtriangle
**********

quietly {
capture program drop rtriangle
program define rtriangle
    syntax [, min(real 0.0) max(real 1.0) mode(real 0.5)]
    local U = runiform()
    if `U' < (`mode'-`min')/(`max'-`min') {
        global x = `min' + sqrt(`U'*(`max'-`min')*(`mode'-`min'))
    }
    else {
        global x = `max' - sqrt((1-`U')*(`max'-`min')*(`max'-`mode'))  
		    }
end
}


***
*pp 
***

quietly {
*The pp program computes the per-person benefit as a function of income.
capture program drop pp
program pp
args Z
*Deal with unacceptable values.
if `Z'<$min {
	global pp=0
}
if `Z'>$max {
	global pp=0
}
*Compute pp for values below the mode.
if `Z'<= $modepp & `Z'>= $min {
	local slope_below=(1-$minpppct)/($modepp-$min)
	global pp=$minpppct+(`Z'-$min)*`slope_below'
}
*Compute pp for values above the mode.
if `Z'<= $max & `Z'> $modepp {
	local slope_above=(1-$maxpppct)/($max-$modepp)
	global pp=1-(`Z'-$modepp)*`slope_above'
}
*di "$pp"
end
}


****
*cdf
****

quietly {
*The cdf program implements the cdf of the truncated triangular distribution (for population) at any value, X. 
capture program drop cdf
program cdf
args X
*Deal with unacceptable values.
if `X'<$min {
	global cdf=0
	*di "Below min"
}
if `X'>$max {
	global cdf=1
	*di "Above max"
}
***triangular
if $unif==0 {
*Compute cdf for values below the mode.
if `X'<= $modepop & `X'>= $min {
	global cdf = (`X'-$min)^2/(($max-$min)*($modepop-$min))
} // close if below mode
*Compute cdf for values above the mode.
if `X'>$modepop & `X'<=$max {
	global cdf = 1-($max-`X')^2/(($max-$min)*($max-$modepop))
} //close if above mode.
} //close triangular
***uniform
if $unif==1 {
	global cdf = (`X'-$min)/($max-$min)
} //close if uniform	
end
}


*******
*simpop
*******

quietly {
*The simpop program computes the normalization factor that accounts for the fact that you don't know maxpp and possibly also total population.
*The first argument is bandwidth, which it gets from the wnb program. The second is the proportion of the population in the band, if you are using the program to compute simpop manually when you have multiple bands.
capture program drop simpop
program simpop
args bandwidth nj
*Set nj to 1 if computing for a single band (in which case nj will not be specified)
capture local t=`nj'
if _rc!=0 local nj=1
*Set up macro that tallies the simulated unweighted net benefit
global nbprop = 0
*Loop over bands.
forvalues i = 1/$k  {
*Define macros for the lower and upper bounds of the band.
	local lb = $min+(`i'-1)*$bw
	local ub = `lb'+$bw
*Compute the cdf at the upper and lower bounds and the mass between.
	cdf `lb'	
	local cdf_lb=$cdf
	cdf `ub'
	local cdf_ub=$cdf
	local mass = `cdf_ub'-`cdf_lb'	
*Compute pp at the midpoint of the band
	local temp=(`lb'+`ub')/2
	pp `temp'
*Compute the net benefit in proportional terms
	if "$pp_units"=="0" global avgwtp_simpop=1
	local nbpropi = `mass'*$pp*`nj'*$avgwtp_simpop
*Tally total net benefit in proportional terms
	global nbprop = $nbprop+`nbpropi'
} //close i
*Compute simpop for use in wnb
global simpop = $nb/$nbprop

*Display nbpropj for use when computing simpop for a set of bands with known proportions.
if _rc==0 di "nbpropj = $nbprop
end
}



****
*wnb
****

quietly {
*The wnb program computes the utility-weighted net benefit across a single band. 
*The first argument is an indicator for whether pp is in units. If it is 
*The second argument is simpop. If it is left blank, the program functions as if the band is the entire population, computing simpop using the simpop program. If the band is one of a set of bands, use the manually computed simpop as the second argument.
capture program drop wnb
program wnb
args pp_units simpop
noisily di "here3"
global pp_units=`pp_units'
capture local t=`simpop'
if _rc==0 global simpop=`simpop'
*Loop over different numbers of sub-bands. 
forvalues k = $bnmin(2)$bnmax {
global k=`k'
*Define sub-band width (bw)
global bw = ($max-$min)/`k'
*Compute simpop if it is not given. (I.E. if you have a single band rather than one of a set of bands.)
simpop $bw //This is supposed to be triggered by if _rc!=0 but that wasn't working.
*noisily di "simpop=$simpop"
*Set up macro that tallies the weighted net benefit
global wnb = 0
*generate eqfac as a function of epsinc
preserve
use "C:\Users\aclan\Dropbox\x Work\2 Research\Distributional concerns\Practicalities\Equivalization for MC.dta",clear
capture drop temp
gen temp=equiv_hh^$epsinc
save "C:\Users\aclan\Dropbox\x Work\2 Research\Distributional concerns\Practicalities\Equivalization for MC.dta",replace
restore
*Loop over sub-bands.
forvalues i = 1/`k' {
*Define macros for the lower and upper bounds of the sub-band.
	local lb = $min+(`i'-1)*$bw
	local ub = `lb'+$bw
*Compute the cdf at the upper and lower bounds and the mass between.
	cdf `ub'
	local cdf_ub=$cdf
	cdf `lb'	
	local cdf_lb=$cdf
	local mass = `cdf_ub'-`cdf_lb'
	*di "mass=`mass'"
*compute ppprop_i
	pp (`lb'+`ub')/2
	local ppprop_i=$pp
	*di "ppprop=`ppprop_i'"
*Compute wtp for the sub-band
	if "`pp_units'"=="0" local wtp=1
	else {
	local wtp_lb=$avgwtp*(`lb'/$avgwtpinc)^$epswtp
	local wtp_ub=$avgwtp*(`ub'/$avgwtpinc)^$epswtp
	local wtp=(`wtp_lb'+`wtp_ub')/2
	} // close else
*Compute nbprop_i x simpop, possibly multipled by wtp
	local impact = `mass' * `ppprop_i' * $simpop * `wtp'
*Compute the weights 
	global lb=`lb'*$alpha/`lb'^$beta
	global ub=`ub'*$alpha/`ub'^$beta
	local mp=(`lb'+`ub')/2
	global mp=`mp'*$alpha/`lb'^$beta
	preserve 
	use "C:\Users\aclan\Dropbox\x Work\2 Research\Distributional concerns\Practicalities\Equivalization for MC.dta",clear
	quietly sum temp [weight=hhwt] if inrange(hhincome,$lb-1000,$lb+1000) $hhdist
	local fac=r(mean)
	global wgt_lb=($med/$lb)^$epsinc*`fac'
	quietly sum temp [weight=hhwt] if inrange(hhincome,$ub-1000,$ub+1000) $hhdist
	local fac=r(mean)
	global wgt_ub=($med/$ub)^$epsinc*`fac'
	quietly sum temp [weight=hhwt] if inrange(hhincome,$mp-1000,$mp+1000) $hhdist
	local fac=r(mean)
	global wgt_mp=($med/$lb)^$epsinc*`fac'
	restore
	local wgt=(($wgt_lb+$wgt_ub)/2+$wgt_mp)/2
	if `wgt'>$thresh_hi local wgt=$thresh_hi
	if `wgt'<$thresh_lo local wgt=$thresh_lo
*Compute the weighted impact for the band and add it to the cumulative tally.
	local wgtd_impact=`impact'*`wgt'
	global wnb = $wnb+`wgtd_impact'
} //close i
} // close k
end
}

}